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One major task in clinical pharmacology is to determine the pharmacokinetic-pharmacodynamic 
( PK-PD ) parameters of a drug in a patient population. NONMEM is a program commonly used 
to build population PK-PD models, that is, models that characterize the relationship between a 
patient 's PK-PD parameters and other patient specific covariates such as the patient 's (pat ho) phy- 
siological condition, concomitant drug therapy, etc. This paper extends a previously described 
approach to efficiently find the relationships between the PK-PD parameters and covariates. In a 
first step, individual estimates of the PK-PD parameters are obtained as empirical Bayes estimates, 
based on a prior NONMEM fit using no covariates. In a second step, the individual PK-PD 
parameter estimates are regressed on the covariates using a generalized additive model. In a third 
and final step, NONMEM is used to optimize and finalize the population model. Four real-data 
examples are used to demonstrate the effectiveness of the approach. The examples show that the 
generalized additive model for the individual parameter estimates is a good initial guess for the 
NONMEM population model. In all four examples, the approach successfully selects the most 
important covariates and their functional representation. Tlte great advantage of this approach is 
speed. The time required to derive a population model is markedly reduced because the number 
of necessary NONMEM runs is reduced. Furthermore, the approach provides a nice graphical 
representation of the relationships between the PK-PD parameters and covariates. 

KEY WORDS: pharmacokinetics; population analysis; model building; generalized additive 
models; NONMEM. 
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INTRODUCTION 

Studies of population pharmacokinetics and pharmacodynamics have 
led to the appreciation of the great degree of variability in pharmacokinetic- 
pharmacodynamic (PK-PD) parameters across patients. Numerous studies 
have quantified the effects of factors, such as age, gender, disease states, 
concomitant drug therapy, etc., on the pharmacokinetics and pharmacodyn- 
amics of drugs, with the goal of accounting for interindividual variability. 
NONMEM is a widely used program for pharmacokinetic and pharmacody- 
namic population analysis (1). NONMEM describes interindividual vari- 
ability in terms of fixed and random effects. The fixed effects relate PK-PD 
parameters to covariates. The interindividual random effects quantify the 
residual unexplained variability. An additional random effect quantifies 
intraindividual and measurement variability. 

Finding a population model that adequately describes the data can be 
a painstaking and complicated task. This task comprises not only finding 
the covariates that significantly influence the PK-PD parameters but also 
determining the shape of the relationships between covariates and param- 
eters. In particular when the PK-PD models become complicated and the 
effects of numerous covariates must be considered, the number of possible 
models increases dramatically. Obviously, one cannot try all possible models 
because the computer time required would be unacceptable. 

In this paper we describe an extension of a previously suggested 
approach to efficiently find the relationships between covariates and PK-PD 
parameters (2). Four real-data examples are used to illustrate the stepwise 
procedure to building a population model that adquately describes the data. 
In describing the approach, we have divided the paper in several sections. 
The general procedure of the analysis is described in the first section. In the 
Data section, the real-data examples are outlined; the Application section 
gives the application of the analysis to the real-data examples, followed by 
a short Discussion section. 



GENERAL PROCEDURE OF NEW APPROACH 

The data analysis is performed using a stepwise approach. Briefly, indi- 
vidual empirical Bayes estimates of PK parameters are obtained. Subse- 
quently, the relationship between the individual PK parameter estimates and 
covariates is modeled using a generalized additive model (GAM) to allow 
nonlinear covariate-parameter relationships to be discovered. Finally, the 
population model is built using NONMEM on the basis of the GAM analy- 
sis of the previous step. The procedure is now described in greater detail. 
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1. In a first step, a basic population PK model without any covariates 
is estimated. In all examples the NONMEM (Version 4, level 1.0) model for 
the yth concentration measurement in the iih individual is given by 

Q=/(Py, '</)(!+ *,) (i) 
where / is the model predicting the jth concentration (usually a mono- or 
biexponential) and p, = (/>■,, p 2i , • * . , p,m) are the m PK parameters of the 
/th individual. £, y represents the residual departure of the model from the /th 
observation available from individual /. The kth element of p, is modeled as 

p ki =8 k exp(r] ki ) (2) 

where 0 k is the population mean of/?*, and exp (tj ki ) expresses the (random) 
difference between Q k and p ki . { e 0 } are assumed independent normally distri- 
buted with mean zero and variance a 2 . 7] 2i , . . . , ij m ,-) are assumed 
to be independent multivariate normally distributed, with mean zero and 
with variance-covariance the diagonal matrix Q with diagonal elements 
(al . . . , (ol). We assume no covariance between elements of rj at this stage 
to be sure that each covariate has the opportunity to appear to be related 
to each r/. Subsequently, empirical Bayes estimates of the p ki are obtained 
by standard methods using the estimated values of 0 k , a 2 and Q. (3,4). 

2. With the estimates of p ki from Step 1 treated as "data," this step 
corresponds to the now classical regression problem of variable selection, 
and we take advantage of the considerable recent work done on this problem 
by others (e.g., 5). In this step, a regression model is derived to model the 
dependence of the individual PK parameter estimates p ki on the covariates 

X\ i j %2i » • • • > Xni • 

^=^w,...,X m ) (3) 

where g k is a multidimensional smooth function. Without some constraints, 
such a genera] description of the data would entail serious problems. First, 
there is the problem of interpreting the function in dimensions higher than 
2. Second, there is the problem, as Bellman (6) put it, of the "curse of 
dimensionality" : The number of points required to fill a space to a given 
density grows exponentially with the dimension of the space but usually the 
data do not. This makes the estimation of complex multivariate functions 
using sparse noisy data highly inaccurate. An often-used simplification is 
to characterize the relationship between parameters and covariates using a 
multiple linear regression model 

n 

Pki=ako+ £ a kl X ti (4) 
/-i 

However, this model makes a strong assumption about the dependence of 
p k on Xu namely, that the dependence is linear in each of the predictors. 
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A more general recent approach is the use of so-called generalized additive 
models (GAM) (5). 

/>*, = <*«,+ £ gk&Xi-i (5) 

where gkt(X/,) is an arbitrary univariate function with £^ , gt/CX//) = 0. The 
functions gkt(Xa) can be represented by any function. The use of spline 
functions (for description of splines see Appendix) is very appealing, because 
they are flexible, yet easy to use. In a subsequent step the spline functions 
can be replaced by a parametric representation, if appropriate. By assuming 
an additive structure, GAM allow straightforward interpretation and display 
of the role of the various covariates. The price one pays for this simplification 
is that interactions between covariates are not included. Nonetheless, a GAM 
is often a useful approximation to the necessarily more complex true multi- 
variate regression surface. 

The GAM is built using a stepwise addition/deletion method (7). This 
method steps through a series of models along a path determined as follows: 
Each of the covariates is allowed to be left out of the model at each step, or 
to enter the model in any of several prespecified functional representations. 
At each step the model is changed by addition (deletion) of the single term 
that results in the largest decrease in the Akaike information criterion (AIC) 
(8). In this context, the AIC is proportional to the residual sum of squares 
from the GAM fit, but adds a penalty, proportional to the number of param- 
eters in the model. The search is stopped when the AIC has reached a 
minimum value. Calculations were performed using the statistical program 
S-plus (Version 3.0, Statistical Sciences Inc., Seattle, WA). A Fortran pro- 
gram for GAM fits is also available elsewhere (ref. 5,- p. 307). This step 
represents our elaboration on the method of Maitre et al. (2) who simply 
advise examination of graphical displays at this step so as to proceed to the 
next and last step. 

3. In a third and final step a parametric or semiparametric mixed effects 
model describing the relationship between covariates and PK parameters is 
built using NONMEM. The regression models found in the second step 
serve as an initial guess for the final population model. In the NONMEM 
analysis the addition of covariance terms between the PK parameters is 
explored. Plots of the individual PK parameter estimates from Step 1 versus 
each other serve to indicate which covariance terms should be included. 

DATA 

Data on four drugs are discussed. They are disguised for confidentiality. 
The analyses reported here are not meant to be definitive but are presented 
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merely to illustrate the use of the model-building method described above. 

Drug A is a nonsteroidal anti-inflammatory that was administered as a 
single oral dose of 5 mg/kg or 10 mg/kg to 93 children with febrile illness. 
Two to six plasma samples were obtained from each individual up to 8 hr 
postdosing (total of 41 1 plasma samples). The following demographic data 
were recorded: gender (SEX), male (53), female (40); race (RACE), 
Caucasian (17), black (76); location (LOC), clinic (65), hospital/inpatient 
(28); food since dosing (FED), no (25), yes (68); dose level (DRG), 5 mg/ 
kg (45), 10 mg/kg (48); height (HT), 64-148 cm; weight ( WT), 5.6-45.5 kg; 
age (AGE), 3-132 months. 

Drug B is an antiarrhythmic, administered orally to 136 hospitalized 
men for various arrhythmias. Plasma samples were obtained for routine 
clinical purposes (total of 361 samples). The following demographic data 
were recorded: race (RACE), Caucasian (91), Latin (35), black (10); 
smoking (SMO), no (91), yes (45); ethanol abuse (ET), no ethanol abuse 
or social drinker (90), current ethanol abuse (16), history of ethanol abuse 
(30); congestive heart failure (HF), no or mild (56), moderate (40), severe 
(40) ; creatinine clearance (RF), < 50 ml/min (48), >50 ml/min (88) ; weight 
(WT), 41-119 kg; height (HT), 154-202cm; age (AGE), 42-92 years; 
a- 1 -acid glycoprotein concentration (GLP), 39-316 mg/dl. 

Drug C is an antihypertensive, administered twice daily in a dose range 
of 1 mg/day up to 20 mg/day in patients with mild, moderate, or severe 
hypertension. Each patient started with the lowest dose and the dose was 
increased until the patient's arterial pressure was controlled. After at least a 
period of 8 weeks continuous therapy the pharmacokinetics of Drug C were 
studied in 64 patients during a 12-hr dosing interval. After 4 weeks, 35 
patients were restudied. A total number of 887 plasma samples were 
obtained. The following demographic data were recorded: gender (SEX), 
male (63), female (36); race (RACE), Caucasian (67), black (32); visit 
number ( VIS), visit 15 (53), visit 16 (46) ; Smoking (TOE), no (74), yes (25) ; 
prior therapy (PT), no (3), yes (96); cotreatment with hydrochlorothiazide 
(HCTZ), no (44), yes (55) ; cotreatment with propranolol (PROP), no (85), 
yes (14); other concomitant therapy (CON), no (24), yes (75); food (FF), 
fed (50), fasting (49); height (HT), 140-188 cm; weight (WT), 51-139 kg; 
age (AGE), 24-69 years; serum creatinine (SECR), 0.6-1.8 mg/dl. 

Drug D is a broad spectrum antibiotic. The pharmacokinetics of drug 
D were studied in 74 critically ill patients with infections with organisms 
sensitive to drug D. The drug was administered twice daily, 200 mg or 
400 mg, in 1-hr infusions. For each individual three plasma samples were 
obtained during a 12-hr dosing interval. In total 113 dosing intervals were 
studied (337 plasma concentrations). The following demographic data were 
recorded: age (AGE), 18-84 years; weight (WT), 43-125 kg; creatinine 
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clearance (CLCR), 0.4-312 mL/min; Glasgow score (GLAS), 3-15; simpli- 
fied acute physiology score (SAPS), 1-26; albumin (ALB), 17-40 g/L; bili- 
rubin (BIL), 4-150 jimol/L; alanine amino transferase (ALAT), 3- 
200IU/L; alkaline phosphatase (AP), 32-615 IU/L; prothrombin level 
(PT), 36-100% normal; systolic blood pressure (SPB), 60-175 mmHg; heart 
rate (HR), 60-170 beats/min ; artificial ventilation (A V), no (42), yes (71); 
gender (SEX), male (86), female (27); center where study was performed 
(CEN), center 1 (44), center 2 (69); administration of high dose catechol- 
amine (CAT), no (92), yes (21). 

APPLICATION 

The stepwise procedure of the new approach is demonstrated in full 
using drug A. The other examples are used to address certain specific issues. 
The purpose of the first step is to determine the basic PK model and to 
obtain empirical Bayes estimates of the model parameters for each individual 
using NONMEM. A simple one-compartment model with first-order 
absorption best described the data of Drug A. 

Changing Covariates Within One Individual 

A complication may arise with the next step, the empirical Bayes estima- 
tion of the model parameters, when the value of a covariate changes in time 
within one individual. Actually, what one wants is an estimate of the model 
parameters for each constant value of a covariate within one individual. 
However, some covariates may change literally from one measurement to 
the next, and one cannot obtain precise empirical Bayes estimates of the 
model parameters for each separate measurement. To solve this problem 
each subject's data is subdivided into disjoint contiguous time periods in 
each of which the covariates are fairly constant. The data corresponding to 
each of these time periods are considered to be a new "individual." The 
length of these time periods depends on the number of samples taken, the 
variability of the covariates, and complexity of the PK model. In this way 
it is possible to obtain empirical Bayes estimates of period-specific model 
parameters which can be correlated to the corresponding time averaged 
values of the covariates during the period. Sometimes this subdivision comes 
naturally, e.g., when the subjects are studied during different dosing intervals 
and several samples are taken during each dosing interval. This is the case 
for Drugs C and D. For those drugs, in fact, some of the covariates changed 
markedly between the dosing intervals studied, but were constant within the 
short period of the dosing interval. For each dosing interval, on average, 9 
(Drug C) and 3 (Drug D) concentration measurements were available, 
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allowing a reasonable empirical Bayes estimate of the individual parameters. 
For Drug B one of the covariates (GLP) changed markedly with each con- 
centration measurement for each individual. To proceed with these data 
individual records were not subdivided by period because on average only 
2.7 samples were available per individual. Instead, the time-averaged value 
of GLP (and other covariates) was used in Step 2. 

The 0, Q, and cr required for the empirical Bayes estimation were 
obtained by fitting the basic PK model to the modified data sets (with 
subdivisions when necessary) for drugs A-D. At this stage, parameters p ki 
and pa, l^k were assumed to be uncorrected, as previously used in the 
definition of Q. 

Initial Screening for Covariate Effects 

Before fitting a GAM by stepwise addition/deletion, the individual PK 
parameter estimates are regressed independently on each covariate, which is 
equivalent to the method of Maitre et ai (2). This initial screening gives a 
first impression of the relative importance of several covariates (i.e., their 
ability to reduce the residual sum of squares) and of the shape of the relation- 
ships between covariates and PK parameters. The results of the initial screen- 
ing for drug A are summarized in Table I, using the parameter clearance 
(CI) as an example. The covariates were modeled linearly or as a natural 
cubic spline with two internal knots at the 33 and 66% quantiles (equivalent 
to 4 df ). The significance (p) of the linear models was tested against the 
constant model using the Ftest. The RSS of the constant model was 12.15. 
The significance of the splines was tested against the linear models. For drug 
A the initial screening resolved that (/7<0.05) AGE, WT, HT, DRG, and 
FED significantly influenced CI. The relationships between CI and the con- 
tinuous covariates (AGE, WT, HT) were best described by a linear model. 



Table I. Change in Residual Sum of Squares (ARSS) 
After Independent Correlation of Each Covariate with 
the Clearance of Drug A 



Covariate 


Linear vs. constant 


Spline vs. linear 


&RSS 


P 


ARSS p 


WT 


-3.89 


<0.01 


-0.10 0.57 


AGE 


-3.88 


< 0.01 


-0.13 0.49 


HT 


-3.64 


< 0.01 


-0.21 0.33 


FED 


-0.61 


0.03 




DRG 


-0.58 


0.03 




RACE 


-0.46 


0.06 




LOC 


-0.07 


0.45 




SEX 


-0.03 


0.60 
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The use of a natural cubic spline to describe these relationships did not result 
in a significant improvement of the fit. 

GAM Analysis 

The GAM is then built using stepwise addition/deletion. At each sub- 
step of the GAM steps, each of the covariates is allowed to be left out of 
the model, to enter linearly, or to enter as a natural cubic spline with two 
internal knots at the 33 and 66% quantiles. For drug A the GAM analysis 
indicated that CI is a linear function of WT and a function of the categorical 
covariates DRG and RACE, and Kis a linear function of WTand a function 
of the categorical covariates LOC and DRG. Table II summarizes the path 
taken to these models and shows some of the other GAMs that are close to 
the best fit one. In the notation of the model equations, ~ means "is a 
function of" and 4- signs are to be interpreted figuratively and not literally. 
In each step the term that results in the largest decrease in the AIC is added 
to the model. The residual sum of squares (RSS) and change in RSS from 
the previous model are also shown. WT has the greatest influence on CI, 
and it alone explains 32% of the variability in C/, while all three factors 
( WT, DRG, and RACE) explain 38% of the variability. HT and AGE, which 
appeared to be important factors in the initial, one covariate at a time, 
screening (/?<0.01), do not appear in the additive model. The reason for 
this is that these two covariates are highly correlated with WT (r> 0.940) 



Table II. Path Taken to the Final GAM for the Clearance (CI) and Volume (V) Parameters 

of Drug A 







CI model 0 










V model 0 






Step 


Term 


&RSS 


RSS 


AIC 


Step 


Term 


ARSS 


RSS 


AIC 


1 


constant 




12.16 


234.3 


1 


constant 




40.1 


345.3 


2 


+ WT 


-3.89 


8.26 


200.5 


2 


+ HT 


-6.6 


33.6 


330.7 


3 


+ DRG 


-0.36 


7.90 


198.2 


3 


+ LOC 


-4.2 


29.3 


320.3 


4 


+ RACE 


-0.32 


7.58 


196.3 


4 


+ DRG 


-1.1 


28.2 


318.6 



Summary of models close to the minimum model 
CI models V models 





Equation 6 


AIC 


Equation 6 


AIC 




WT+ DRG + RACE 


196.3 


V-WT+LOG + DRG 


318.5 


c/~ 


WT + DRG +RACE+ LOC 


197.3 


V~HT+LOC+DRG 


318.6 


c/~ 


WT+ DRG + RACE+ SEX 


197.3 


V~ WT+ LOC+ DRG + DRG + SEX 


319.7 








V~ HT+ LOC + DRG + SEX 


319.9 








V~ HT+ LOC+ DRG + FED 


320.0 



"In each step the term that results in the largest decrease in the AIC is added (+) to the model. 
*In the notation of the model equations ~ means "is a function of* and the + signs are to be 
interpreted figuratively and not literally. 
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V model 0 


:p Term 


&RSS 


RSS 


AIC 


constant 
+ OT 
+ LOC 
+ DRG 


-6.6 
-4.2 
-1.1 


40.1 
33.6 
29.3 
28.2 


345.3 
330.7 
320.3 
318.6 


ac minimum model 








V models 






Equation* 




AfC 


WT+LOG+DRG 318.5 
HT+ LOC + DRG 318.6 
WT+LOC+DRG+DRG+SEX 319.7 
//7> LOC+ DRG + SEX 3 1 9.9 
HT+ LOC+ DRG + FED 320.0 
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and therefore add nothing to it. Figure I shows a plot of all covariates found 
in the GAM for C/, vs. a variable indicating their effect on CI (see legend 
to Fig.). The spline for WT is displayed instead of the linear relationship 
found in the minimum model. Figure 1 clearly shows that the relationship 
between CI and WT is predominantly linear. The nonlinearity is apparently 
caused by only one subject on the far right of the WT scale. 

During the model search numerous models are visited. The GAM was 
derived using both a forward search starting with the NULL model and a 
backward search starting with the FULL model, which includes all covari- 
ates. In all, 71 different models were tested for CI. For CI the forward and 
backward search resulted in the same model. For V a small difference was 
found: HT appeared in the forward search, while appeared in the back- 
ward search. However, these two covariates are highly correlated which 
makes this understandable. The covariates HT (or WT) and LOC have the 
biggest effect on Kand explain 27% of the variability. The addition of DRG 
results in a minor additional drop in AIC. 

Final NONMEM Analysis 

In the third and final step the population model is built using NON- 
MEM on the basis of the GAM results of the previous step. At a first 
substep, only linear relationships between covariates and parameters are 
considered. In doing so, basically two approaches can be used: (i) Include 
in the initial NONMEM model all covariates that appear in the GAM with 
the lowest AIC value, and in other models close to it, or (ii) start with the 
minimum AIC model as the initial NONMEM model and test later whether 
the other covariates that appear in (Step 2) models close to the minimum 
AIC model significantly reduce —2 log likelihood (—ILL) obtained by the 
NONMEM analysis. (The difference in — ILL is asymptotically x 2 distri- 
buted). The potential additional covariates can be tested according to their 
order of importance in Step 2. For Drug A approach (i) was used, because 
the total number of covariates was relatively small. The following models 
for CI and V were used in the initial NONMEM model (the covariates are 
centered at their medians) 

C/= 0, + 0 2 ( WT- 15) + 6sDRG+0 4 RACE+ 0 S LOC+ 0 6 SEX (6) 

K= 0 7 + 0%(WT- 15) + 0 9 LOC + 0 w DRG+0 u SEX+ 0 X2 FED (7) 

A covariance term between CI and V was added to the model, because the 
plot of the empirical Bayes CI and V estimates showed a strong correlation. 
The parameter estimates found in Step 2 were used as initial estimates for 
the NONMEM fit. In general, a good agreement was found between the 
parameter estimates in the GAM and NONMEM analyses. Including the 
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Fig. 1. Fit of the GAM for the CI of drug A. Each plot is the contribution of a term to the 
additive predictor of CI. The ordinate label is the expression used to specify the contribution 
of the covariate to the model formula in the S-plus language. The ordinate value can be 
thought of as the (partial) effect of the covariate on CI. Each curve has been centered to 
have an average of zero. The points plotted are the partial residuals: That is, the individual 
CI values minus the prediction based on all other covariates. The same scale is used for the 
y axis in all plots so that the relative importance of the covariates can be compared. The 
two dose level (DRG) groups were labeled 1 = 5 mg/kg and 2 = 10 mg/kg, the race (RACE) 
groups were 1 = Caucasian, 2 = black. For these dichotomous covariates, points are randomly 
shifted by small amounts on the abscissa to display them better. 
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various covariates in the model results in a highly significant drop of 238 in 
—ILL when compared to the basic (Step 1) PK model. The initial [Eqs. (6) 
and (7)] model was further refined by deleting certain parameters from the 
model according to the following criteria: (i) Only a small change in —ILL 
when the parameter is left out, (ii) ±1 SE of the parameter estimate includes 
zero, and (iii) the magnitude of the covariate coefficient is small (i.e., the 
covariate has little effect on the PK parameter). This resulted in the following 
models for CI and V in the final NONMEM model 

C/= 0 , + 0 2 ( WT- 1 5) + 0 3 DRG + 0 4 RA CE (8) 
V= 0 5 + 6 6 ( WT- 15) + 0?LOC+ 0*DRG+ 0 9 SEX (9) 

The influence of LOC and SEX on CI and the influence of FED on V 
was not significant. The 0s of these covariates were poorly determined with 
a standard error larger than the parameter estimate. Fixing the 0s of these 
covariates to zero resulted in an increase of only 3.4 in —ILL. Thus, the 
parameters that could be deleted from the initial model appeared to be the 
ones that were the least important according to the covariate selection in 
Step 2. Table III summarizes the change in -2LL when each of the 0s of 
the covariates appearing in the reduced (and final) NONMEM model is set 
to zero, showing the significance of each of the parameters (a change > 6.6 
in — ILL is significant with /><0.01). For C/ f WT seems to be the most 
important covariate in the final population model. This is in accordance 
with the GAM analysis: WT resulted in a large decrease in AIQ while 
addition of the other two covariates resulted in a relatively small additional 
decrease. For V 9 WT and LOC appear to have the greatest effect, which is 
also in accordance with the results found in the GAM analysis. For V an 
additional covariate, SEX, that appeared in a GAM close to the minimum 
model, was significant. Of the four covariates not found in the minimum 
AIC GAM but in models close to it (Table II), only SEX significantly 
influenced the fit. Including the final set of covariates in the NONMEM 
model resulted in a decrease of the interindividual variability from 46 to 



Table III. Change in -ILL When Each of the 0s of the Covariates 
Appearing in the Final NONMEM Model for Drug A is Set to Zero 



CI model 


V model 


Covariate -2LLL (vs. 0 = 0) 


Covariate -2ALL (vs. 0=0) 


WT 115 


WT 77.5 


DRG 16.0 


LOC 50.2 


RACE 9.9 


DRG 8.3 




SEX 14.4 
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29% for CI and from 46 to 28% for V when compared to the basic PK model 
without covariates. 

Validation of Final NONMEM Model 

As a first check on the NONMEM analysis, the relationship between 
WT and CI and WT and V was modeled by a natural cubic spline (see 
Appendix) with two internal knots at the 33 and 66% quantiles (4 df ). This 
nonparametric representation did not result in a significant change in -ILL, 
confirming the results found in the GAM analysis that these relationships 
are predominantly linear. 

Finally, to check whether our final model was really the minimum 
model, all other covariates were included in the model one by one, but none 
significantly influenced -ILL. Another, quicker approach to test whether 
the final model is the minimum model may be the following. Calculate the 
exp (T] kJ ) (the unexplained residual interindividual variability) as empirical 
Bayes estimates using the final NONMEM model Then, regress the 
exp (T] kJ )s on the covariates that are not in the final NONMEM model. This 
was done for drug A and, as expected, none of the remaining covariates 
correlated significantly with the final model-based empirical Bayes estimates. 
We propose this as a final step in the analysis to confirm whether an appro- 
priate final model has been reached. 

Comparison of GAM and NONMEM Analyses 

The GAM analysis apparently provided both the important covariates 
for the final NONMEM model and the shapes of the relationships between 
covariates and PK parameters for example A. An important question regard- 
ing the suggested approach, however, is whether the GAM analysis always 
suggests all the relevant covariates. Although we cannot know that for sure, 
experience with drugs B-D suggests that it is so. Table IV compares the 
covariates that were selected by the GAM approach and the covariates that 
were significant in the population model derived using NONMEM. For all 
drugs the GAM with the lowest AIC included all the important covariates 
or they were included in GAM models close to the minimum model. This 
latter case happened only once: Drug B, covariate RF. Table IV shows the 
results only for the CI parameter of the four drugs, but the same findings 
held for V. The GAM fit also identified the relative importance of the covari- 
ates. In general, the covariates that resulted in the biggest decrease of the 
AIC in the GAM step were also the most important covariates in the NON- 
MEM model (judged by change in -2LL). Conversely, the covariates that 
were of borderline significance in the GAM fit were of minor significance in 
the final NONMEM model or could be deleted from the model. 



Mandema, Verotta, and Sheiner 
:n compared to the basic PK model 



analysis, the relationship between 
led by a natural cubic spline (see 
; 33 and 66% quantiles (4 df ). This 
ult in a significant change in -ILL, 
VI analysis that these relationships 

J model was really the minimum 
in the model one by one, but none 
quicker approach to test whether 
my be the following. Calculate the 
individual variability) as empirical 
MEM model. Then, regress the 
n the final NONMEM model. This 
none of the remaining covariates 
lel-based empirical Bayes estimates, 
ilysis to confirm whether an appro- 



mlyses 

ided both the important covariates 
hapes of the relationships between 
e A. An important question regard- 
whether the GAM analysis always 
>ugh we cannot know that for sure, 
.t it is so. Table IV compares the 
4 approach and the covariates that 
derived using NONMEM. For all 
luded all the important covariates 
;lose to the minimum model. This 
covariate RF. Table IV shows the 
four drugs, but the same findings 
e relative importance of the covari- 
lted in the biggest decrease of the 
important covariates in the NON- 
L). Conversely, the covariates that 
M fit were of minor significance in 
deleted from the model. 



Building Population Models 523 

Table IV. Covariates Found by the GAM and NONMEM Analysis to Significantly 

Influence CI 

Analysis Covariates" 

Drug A 



GAM 


wr 


DRG 


RACE 










NONMEM 


WT 


DRG 


RACE 

Drug B 










GAM 


GLP 


WT 


RACE ET 


RF 6 








NONMEM 


GLP 


WT 


RACE ET 
DrugC 








GAM 


HT 


SEX 


RACE TOB 


AGE 


HCTZ 


SECR 


WT 


NONMEM 


HT 




RACE 

Drug D 


AGE 


HCTZ 






GAM 


CLCR 


BIL 


WT AGE 


AP 


CEN 


SEX 


SPB 


NONMEM 


CLCR 


BIL 


WT AGE 




CEN 




SPB 



"The covariates are listed in their order of appearance in the GAM. 
*Found in a GAM model dose to the minimum AIC model. 



Nonlinear Covariate Relationships 

Another important question is whether the GAM step identifies impor- 
tant nonlinearities in the relationships between covariates and PK param- 
eters. In this respect, a nice example is provided by drug C In the GAM 
with the lowest AfC, HT and AGE were two continuous covariates that were 
of importance. Figure 2 shows the contribution of these two covariates to 
the GAM (the effects of the other covariates are not shown). The relation- 
ships were characterized by a natural cubic spline with two internal knots 




140 150 160 170 160 190 30 40 50 60 70 

ht age 



Fig. 2. Plot of the contribution of HT and AGE to the GAM for C/ of drug C (other 
relevant covariates are not shown). The relationships are described by a natural cubic spline 
with two internal knots at the 33 and 66% quantiles. See legend to Fig. 1 for further details. 
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at the 33 and 66% quantiles. The shapes found in the GAM suggest that some 
nonlinear representation may be appropriate for the NONMEM model. 
Subsequently, in the NONMEM fit a natural cubic spline (with two internal 
knots at the 33 and 66% quantiles, constrained to monotonicity) was used 
to describe the relationship between CI and HT and CI and AGE. Figure 3 
shows the shapes of the spline representations found in NONMEM. Note 
the similarity between the splines found in the GAM analysis and those from 
the NONMEM fit. Using the splines in the NONMEM fit resulted in a drop 
of 30.0 points in — ILL (4 additional parameters) when compared to the fit 
using a linear relationship. The particular shape of the splines suggest that 
HT only has an effect on CI when the subjects are taller than 160 cm. Age 
only appears to influence CI when the subjects are older than 60. Therefore 
the natural cubic splines were replaced by simpler (linear spline) models: CI 
was assumed to be linearly related to HT for subjects taller than 160, and 
constant for smaller subjects; CI was linearly related to AGE for subjects 
older than 60 years and constant for younger subjects. This simplified model 
resulted in about the same -ILL as the original spline fit, but used fewer 
parameters. The nonlinearities in HT and AGE were both statistically sig- 
nificant. Replacing the nonlinear model for AGE by a linear model resulted 
in an increase of -ILL of 18.8. Replacing the model for HT by a linear 
model increased the —ILL by 7.4 points. 

For drug D the GAM analysis selected a nonlinear relationship between 
CI and WT, A similar relationship was found by NONMEM, resulting in a 




140 150 160 170 180 30 40 50 60 70 



ht age 

Fig. 3. Spline (natural cubic spline with two internal knots at the 33 and 66% quantiles, 
constrained to monotonicity) representation of the relationship between CI and HT, and 
CI and AGE for drug C found in NONMEM fit. 
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significantly better fit than the linear model. The correspondence between 
GAM and NONMEM analyses was not always perfect, however. For drug 
B, a nonlinear relationship between CI and GLP was preferred in the final 
NONMEM model but was not detected by the GAM analysis. This was 
probably caused by the previously noted fact that the GLP changed mark- 
edly from one concentration measurement to the next within individuals, 
but because too few data were available to provide an empirical Bayes 
estimate of CI for each value of GLP, the time averaged value of GLP for 
each individual was correlated with the individual CI estimates in the GAM 
step. This undoubtedly obscured the predominantly intraindividual correla- 
tion between the two, which NONMEM eventually discovered. 



DISCUSSION 

NONMEM is a widely used program to determine the PK and PD 
parameters of a drug in a patient population. One of the advantages of 
developing a population PK-PD model is that the influence of certain factors, 
such as (patho)physiological conditions, concomitant therapy, etc., on the 
drug's pharmacokinetics and pharmacodynamics can be included in the 
population model. This allows identification of possibly causal effects that 
can later be confirmed, the design of patient-specific dosage regimens, and 
optimization of drug dosage adjustment using Bayesian approaches. 

Finding a population model that adequately describes the data can be 
a complicated task. Recently, Maitre et al. (2) proposed a method that forms 
the basis of the method presented here. They proposed to look at scatterplots 
of empirical Bayes estimates of the model parameters versus the covariates 
in order to detect which covariates to use in the population model. The 
problem with this approach is that it is valid only if the covariates act 
independently on the PK-PD parameters. In this paper we elaborate on 
Maitre et a/.'s (2) approach and propose a more formal analysis of the 
empirical Bayes estimates in which we regress the individual empirical Bayes 
estimates on the candidate covariates using a generalized additive model 
(GAM). An important feature of the GAM is that it provides a functional 
representation of the relationship between PK-PD parameters and covari- 
ates. Thus using the GAM analysis, a model that can describe the relation- 
ship between PK-PD parameters and the covariates is built outside 
NONMEM. At a later stage, NONMEM is used to optimize and finalize 
the tentative model thus provided. The great advantage of this elaboration 
is speed. The total time required to build the population model is reduced 
dramatically because the number of NONMEM runs is reduced. The elucid- 
ation of the GAM that best describes the relationship between the individual 
estimates of the model parameters and the covariates is done quickly relative 
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to the series of NONMEM runs required to perform the same step. Further- 
more, the data analysis tools used in the GAM step give a nice graphical 
representation of the relationship between model parameters and covariates, 
something heretofore lacking in the NONMEM software. 

The examples presented in this paper have shown that the derived GAM 
is a very good initial model for the NONMEM analysis. Apart from selecting 
the most important covariates the GAM approach also indicates important 
nonlinearities in the relationships between covariates and PK parameters. 
An important feature of the GAM step is building a model that includes the 
contribution of several covariates. This model provides a better starting 
point than the initial, one covariate at a time, screening (which is equivalent 
to the Maitre et al suggestion). The initial, one variable at a time, screening 
tends to select numerous covariates that do not appear to be of importance 
in the final NONMEM model, usually because they carry the same informa- 
tion as one or a combination of the other covariates. What is even more 
important, the initial screening misses some of the covariates that appear to 
be of importance in the final NONMEM model. The reason here may be 
that their effect is overshadowed by some other covariate or that correlation 
among covariates masks the relationship. A nice example is provided by 
drug D. Of the 16 reported covariates, 13 were found to significantly affect 
CI according to the one-at-a-time initial screening (p<0.05, F test). The 
GAM model narrowed this down to 8 covariates of which 6 were finally 
selected by the NONMEM analysis. One of the covariates finally selected 
for the NONMEM model (SPB) was not selected by the initial screening 
but did appear in the GAM. Another example is drug A. The one-at-a-time 
initial screening selected 6 covariates to significantly affect CI (Table I), of 
which 2 appeared in the final NONMEM model (Table III). However, one 
covariate (RACE), which is significant in the final NONMEM model, was 
not selected during the initial screening. Conversely, the GAM analysis 
selected the same 3 covariates as appeared in the final NONMEM model 
(Tables II and III). Similar results were found for the other two drugs. 

A disadvantage of the approach is that its success depends on the quality 
of the individual empirical Bayes parameter estimates. These estimates tend 
to be "shrunken" towards the population mean, especially when few data 
are available from an individual (or period within an individual). This bias 
of the empirical Bayes estimates complicates the detection of relevant covari- 
ates, and perhaps even more, it complicates the elucidation of the shape of 
the relationship between the model parameters and the covariates. However, 
if a covariate is found to be of importance in the GAM analysis, it has a 
greater chance of being clinically relevant. The same probably holds for 
shape effects. If an important nonlinearity exists we surmise that it will either 
show up in the GAM or a linear representation may be an appropriate 



Mandema, Verotta, and Sheiner 



Building Population Models 



527 



d to perform the same step. Further- 
he GAM step give a nice graphical 
en model parameters and covariates, 
3NMEM software, 
er have shown that the derived GAM 
JMEM analysis. Apart from selecting 
A approach also indicates important 
/een covariates and PK parameters, 
is building a model that includes the 
is model provides a better starting 
i time, screening (which is equivalent 
tial, one variable at a time, screening 
it do not appear to be of importance 
because they carry the same informa- 
>ther covariates. What is even more 
ome of the covariates that appear to 
EM model. The reason here may be 
tie other covariate or that correlation 
hip. A nice example is provided by 
13 were found to significantly affect 
ial screening (/7<0.05, F test). The 
3 covariates of which 6 were finally 
hie of the covariates finally' selected 
not selected by the initial screening 
xample is drug A. The one-at-a-time 
} significantly affect CI (Table I), of 
: M model (Table III). However, one 
- in the final NONMEM model, was 
ng. Conversely, the GAM analysis 
;ared in the final NONMEM model 
e found for the other two drugs. 
:hat its success depends on the quality 
neter estimates. These estimates tend 
tion mean, especially when few data 
riod within an individual). This bias 
cates the detection of relevant covari- 
cates the elucidation of the shape of 
imeters and the covariates. However, 
tance in the GAM analysis, it has a 
vant. The same probably holds for 
ity exists we surmise that it will either 
presentation may be an appropriate 



approximation. So far we have observed that the GAM analysis is inclusive. 
For the four drugs a total of 27 covariates were found to be of significance 
in the NONMEM models. Of these covariates 25 were found in the GAM 
with the lowest AIC value for the drugs. The remaining 2 were found in 
models close to the minimum AIC model. Thus, apparently the shrinkage is 
not a major problem in detecting important covariates. As mentioned in the 
Data section, an additional problem is that interactions between covariates 
are not included in the GAM analysis. Recently, two new approaches to 
covariate and shape selection : PI (9) and MARS (10) have become available. 
These approaches automatically search for combinations of covariates show- 
ing important interactions and determine their functional representation 
using splines. We are currently investigating the use of these approaches at 
Stage 2 and the implementation of the derived models at Stage 3. 

A final note of caution : In model building subjective judgment plays a 
much larger role than in other more formal pursuits. The final model is not 
the only, or true one; only a useful one, and although our experience is 
encouraging, further work with the approach we suggest is required to deter- 
mine its true utility. 
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APPENDIX 

A spline function, (sp(x), where jc is a predictor variable), is charac- 
terized by a strictly increasing sequence of real values (b } < • * * < b n ) called 
knots, and by its order k. In the interval (6, , 6,+i) the spline is a polynomial 
of order k - 1 , and across each 6, , derivatives of order 0 to k - 2 are continu- 
ous. All spline calculations in NONMEM were made using the package (B- 
spline) based on PPPACK (11). B-spline is available from D. Verotta, Box 
0626, University of California-San Francisco, San Francisco, CA 94143, 
upon request (email david@c255.ucsf.edu or david@ucsfc255.bitnet). 
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